The Annals of Applied Statistics

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Estimating the correlation between subject-specific rates of change is a common objective in longitudinal studies. A natural approach is to fit ordinary least squares regressions per subject and correlate the resulting slope estimates directly. Under complete data, this slopes-first estimator is algebraically equivalent to a contrast-based plug-in estimator derived from the joint variance-covariance matrix of the two signals across time. Under independent dropout, however, the slopes-first approach implicitly applies heterogeneous weights across missing-data patterns, introducing bias whenever the cross-covariance between the two signals is not constant over time. We formalise this result and propose an unbiased alternative that estimates the joint covariance matrix before applying the contrast. We further extend our results to include covariates as well. Simulations confirm that the slopes-first estimator is biased under non-constant cross-covariance while the proposed alternative remains unbiased. Further, we show in a real-world example that the two estimators can lead to different scientific conclusions and recommend that researchers working with longitudinal data inspect the missing data pattern to make an informed decision about the estimator of choice.

Spatial expansion is a common feature of many ecological processes, including biological invasions, range shifts, population recovery, and the spread of pathogens and human activities. Although such expansion often arises from diffusive processes, there remain few general statistical methods to quantify these dynamics. Here, we present a diffusion-enhanced spatio-temporal model (DESTM) that approximates diffusive process using Gaussian Markov random fields (GMRF). The model works by incorporating the statistical interaction between spatial diffusion and temporal lags into a joint precision structure of the latent field. To demonstrate, we simulate data from a diffusion movement process using a continuous-time Markov chain (CTMC), and then approximate this mechanistic process using the DESTM. We show that separable spatio-temporal models can lead to large bias in estimates of the local density. In contrast, the DESTM did not appear to bias inferences, and reverts to the separable model when diffusion is not apparent. Finally, we apply the model to two large-scale examples of spatial expansion: (1) the historical spread of the Japanese distant-water longline fleet across the Pacific, representing anthropogenic spatial expansion in the ocean, and (2) the recovery and geographic expansion of Bald Eagle populations across North America. Estimates from both case studies, supported by model comparisons and ecological plausibility, suggest that the DESTM provides a better description of the underlying expansion dynamics, although improvements in forecast performance were more apparent for the Japanese longline fleet. We highlight that the model provides a fast and efficient framework for modeling spatial expansion dynamics, with the potential to represent increasingly complex ecological dynamics through future extensions. We propose that future research explore this framework as a quantitative tool to explain how spatial expansion occurs.

Longitudinal non-Gaussian data are common in biomedical and microbiome studies, where responses may not be normally distributed. However, many existing dynamic graphical models are developed for Gaussian responses, while many non-Gaussian graphical models focus on cross-sectional dependence and do not explicitly estimate temporal transition effects. We propose a latent Gaussian dynamic graphical model for longitudinal non-Gaussian responses, including binary, count, and compositional outcomes. The proposed framework jointly estimates the contemporaneous conditional dependence network and the lag-one temporal dependence structure. We further develop a penalized Monte Carlo expectation maximization procedure for estimation. Simulation studies across multiple graph structures and outcome types show that the proposed method can recover both dependence structures, with more stable recovery for the lag-one temporal dependence structure and more structure-dependent recovery for the contemporaneous conditional dependence network. We apply the proposed method to longitudinal data from the Integrative Human Microbiome Project (HMP2) using a multinomial observation model. The results reveal persistent lag-one temporal dependence and contemporaneous conditional dependence among dominant microbial taxa, demonstrating the utility of the proposed framework for dynamic network analysis of longitudinal non-Gaussian biomedical data.

Somatic mutations in sequenced tumors are often detected with imperfect sensitivity. They can be missed because tumor heterogeneity results in some mutations being absent in a sample at hand or present at very low frequencies which, combined with technical limitations like insufficient sequencing depth, reduces their detectability. It becomes apparent with multi-regional sequencing that some mutations confined to specific areas can be overlooked in single-region analyses. Adjusting for such sensitivity becomes critical when mutation frequencies between groups of subjects are compared, such as responders and non-responders to a given therapy. Otherwise, the estimate of the mutation prevalence, or of the difference in prevalences between the groups, can be biased, leading to inflated type I error. Here, we propose a weighted approach to adjust for imperfect sensitivity of the mutation detection and compared it to the unweighted approach. We assume that a surrogate of mutation sensitivity is observed only among patients with called mutation, and we derive its unconditional mean. Such adjustment leads to less biased estimates of mutational frequencies and better control of the type 1 error.

Natural killer (NK) cells display striking variability in their cytotoxic responses, but whether this variation reflects random events, stable differences between cells, or effects of previous encounters remains unclear. We developed a Bayesian framework that uses single-cell interaction histories to disentangle these sources of variability. The framework models target encounters and killing decisions to quantify population heterogeneity, determine the sample sizes needed to distinguish competing mechanisms, and separate stable cell-to-cell differences from history-dependent behaviour. Synthetic data established when these mechanisms can be reliably distinguished in practice. We tested the framework using time-lapse imaging of NK cells exposed to rituximab or CC-96673, an antibody co-targeting CD20 and CD47. Despite similar overall killing, the framework identified distinct underlying responses: rituximab mainly increased mean killing rate, whereas CC-96673 reduced cell-to-cell variation. Event-count analysis supported continuous population heterogeneity, while ordered interaction histories revealed that both stable differences in killing propensity and previous encounters shape cytotoxic decisions. To make the framework directly usable, we introduce BARRACUDA, an open-source web platform and Python package implementing these analyses, including donor-aware extensions, for single-cell cytotoxicity datasets.

Predicting implantation outcomes in vitro fertilization (IVF) is challenging because clinically relevant information becomes available at different points during treatment. In this study, we developed a stage-aware, leakage-conscious, and explainable machine-learning framework to predict implantation using a public IVF dataset. {beta}-HCG was excluded to avoid potential downstream target leakage. The data were divided into training, validation, and independent test sets. Several machine-learning models, including Logistic Regression, Random Forest (RF), Extra Trees, XGBoost, and Hist Gradient Boosting, were optimized using stratified cross-validation. A fixed RF model was then evaluated at three stages of the IVF process, using baseline couple characteristics, gamete and fertilization variables, and embryological variables, respectively. Model performance improved as additional treatment-related information became available, with the largest improvement observed after the inclusion of embryological variables. The final RF model showed good predictive performance on the independent test set. The number of embryos was consistently identified as the most important predictor across RF, permutation, and SHAP analyses. Overall, the findings suggest that IVF implantation prediction becomes more accurate as later-stage treatment information is incorporated, particularly information related to embryological characteristics. The proposed approach also provides interpretable predictions while reducing the risk of information leakage. Further evaluation using independent clinical datasets is needed before the framework can be considered for clinical use.

Dysfunctions in transcription cause severe human diseases. In eukaryotes, transcription initiation by RNA polymerase (Pol) II requires the formation of the pre-initiation complex (PIC) composed of general transcription factors (GTFs), Pol II, and an essential multiprotein coactivator, Mediator. In vitro experiments with GTFs and Pol II have depicted a linear sequence of PIC assembly. However, PIC assembly paths in vivo at the genomic scale and Mediator's role in these paths remain poorly understood. One strategy characterizes ChIP-seq occupancy alterations in Mediator mutants, but analysis of these data still lacks a formalized interpretation framework. In this study, we address this need by developing a quantitative framework that scores candidate assembly paths according to their ability to explain these alterations. This framework combines a model of PIC assembly based on a system of differential equations with parameter optimization based on a penalized criterion. We applied this new framework to datasets from two yeast Mediator mutants, in the Med17 and Med10 subunits. Our results show that the best-supported PIC assembly paths differ from those inferred in vitro. In particular, we show that it is important to distinguish the core TFIIH module from the TFIIK kinase module, and that similar paths starting with Mediator or TBP receive different levels of support. Furthermore, although multiple paths fit the data equally well, they differ markedly for specific gene groups. Taken together, the results show that the path can depend on the gene, and are also compatible with more than one path being operational in vivo for a single gene. This work improves our understanding of PIC assembly paths and transcription initiation, and opens perspectives for future modeling and experimental studies.

Three-dimensional imaging is rapidly increasing the scale of biological morphometric datasets, but establishing automated anatomical correspondence remains a major bottleneck. Conventional point-cloud registration methods transfer landmarks efficiently, yet poor initialization or high morphological variation can produce geometrically plausible but biologically implausible deformations. Here, we developed a statistical-shape-model (SSM)-regularized coherent point drift (CPD) estimator for automated morphometric landmark transfer. It replaces the unconstrained non-rigid CPD update with optimization in SSM coefficient space and regularizes coefficients by inverse population covariance, so deformations along low-variance biological modes are penalized more strongly. We compared pose-matched conventional CPD, SSM-CPD, and SSM-CPD with refinement in leave-one-out atlas experiments on laboratory mouse skulls (N=249), primate molars (N=115), and squamate vertebrae (N=254) under aligned, mildly perturbed, and severely perturbed initialization. We evaluated these formulations in terms of landmark accuracy, failure-tail behavior, runtime, and preservation of morphometric structure. Under severe perturbation, we also evaluated an integrated recovery procedure that optimizes the starting template, recovers global pose, and then applies SSM-CPD with refinement. Severe misalignment produced widespread failures in pose-matched conventional CPD, whereas SSM regularization reduced median error from 13 - 21% to approximately 1 - 3% of the target bounding-box diagonal. This improvement was highly consistent for skulls and molars. In the morphologically complex vertebral dataset, however, low median error concealed a substantial minority of catastrophic failures (28% of vertebrae above a 6% error threshold), which integrated recovery reduced to 4%. Across datasets, automated landmarks preserved the major axes and pairwise relationships of morphometric shape space, and local CPD refinement improved covariance preservation. Embedding learned biological variation directly into registration broadens the range of starting conditions and anatomies for which automated landmark transfer remains reliable. Our results also show that failure tails and downstream morphometric fidelity must be assessed alongside median accuracy when automated landmarking methods are evaluated.

Understanding the causal pathways linking the gut microbiome to downstream biomarkers and clinical outcomes is central to microbiome research. When the microbiome acts as the exposure, the high dimensionality and compositional structure of the data introduce unique methodological challenges for mediation analysis. To address these challenges, we propose a latent variable mediation framework that captures variation in microbial composition through a microbial balance, defined as the log-ratio between two unknown subsets of taxa. This balance serves as a latent scalar exposure and simplifies the estimation of the overall indirect effect at the community level, while simultaneously identifying specific taxa that contribute to the overall direct and indirect effects. We describe the model's estimation and inference, and illustrate the method using data from the Multiethnic Cohort--Adiposity Phenotype Study, where we examine the mediation pathway from the gut microbiome through lipopolysaccharide-binding protein, a marker of metabolic endotoxemia, to percent liver fat.

In this paper, we introduce Nagini3D, a 3D biological image segmentation method that hybridizes two high-performing yet highly complementary segmentation paradigms: deep learning-based approaches and snake-based methods. Our approach leverages both the high efficiency of data-driven techniques and the convenient shape representation provided by active surfaces, which proves highly practical for incorporating a priori knowledge as well as for the post-analysis of segmented objects. By considering a complete 3D volume, our method ensures segmentation consistency across all three dimensions, unlike many approaches that handle volumes in a slice-by-slice manner. Additionally, it theoretically enables the representation of any object with spherical topology, overcoming the limitations of common 3D methods that rely on object convexity assumptions. Furthermore, the continuous representation of objects allows for a resolution-independent description, facilitating the extraction of various geometric features of interest (e.g., local curvature, principal deformations). We demonstrate both the quality of the segmentation produced by our method and the richness of the accessible geometric features through an extensive set of experiments conducted on multiple real and simulated datasets, including objects with challenging geometries.

Fishery-dependent indices of abundance are often used to inform stock assessments, but fishing effort is generally not randomly distributed. Fishers concentrate sampling in locations expected to be profitable, leading to a form of preferential sampling that can bias model-based indices of abundance. We show that when the strength of preferential sampling is constant over time, it biases only the magnitude of an index, which is absorbed by estimated catchability parameters in most stock assessments. A more interesting case arises when the strength of preferential sampling changes over time; in these cases the index trend can be biased in a way that conventional catchability estimation cannot correct. Using simulations structured on long running trawl survey data from the west coast of the USA for shortspine thornyhead (Sebastolobus alascanus), we evaluate two structurally different methods for correcting time-varying trend bias: a joint spatiotemporal model of the density and sampling processes, and a discrete-choice model of fisher location choice. Our simulations include two mechanisms generating preferential biases: effort tracking local density, and effort responding to a spatially structured, but unobserved anomaly that also inflates catch. We show that the ability of the estimation methods to correct bias depends on the mechanism responsible for concentrating effort; while the joint model is robust to both data-generating mechanisms, the discrete choice model is only able to correct biases when the mechanism is anomaly-driven. As a real-world case study, we apply these corrections to simulated preferentially sampled data and include indices in the shortspine thornyhead stock assessment model. We found that uncorrected time-varying preferential sampling positively biased estimates of a major recruitment event and, to a lesser extent, unfished recruitment. This bias propagated through to bias in overfishing limits that grew as the population was projected forward. While bias-corrected indices still produced some positive biases in overfishing limits, these were substantially reduced compared to the naive approach. As the mechanism of preferential sampling is rarely known in advance, our results suggest that the joint model represents a robust default for correcting fishery-dependent indices of abundance.

Forecasting and causal inference in complex, non-linear dynamical systems are fundamentally challenging due to noisy and high-dimensional time-series data. While Empirical Dynamic Modeling addresses this via state-space reconstruction in Euclidean space, it suffers from the "curse of dimensionality," distorting distance measurements as dimensions increase. This study proposes Forest-Weighted S-map (FORWS) and Forest-Weighted Causal Inference (FORWC), utilizing adaptive "forest weights" derived from random forest ensembles instead of Euclidean metrics. Using simulated and empirical ecological datasets, I demonstrate that these methods exhibit comparable or improved forecasting skill relative to conventional tools, while offering substantial resilience against dynamic process noise and mitigating the curse of dimensionality. Notably, FORWC enables multimodal causal inference, successfully extracting directed interactions between high-dimensional acoustic vectors and scalar temperature data monitored in a honeybee hive. By bridging chaotic physics with machine learning, this framework accurately extracts intrinsic manifold topology, providing a robust paradigm for uncovering causal networks in complex real-world environments.

Selection bias is a common concern in causal inference. Two mechanisms can give rise to selection bias: collider bias or stratifying on an effect modifier of the exposure-outcome causal effect. Research into the impact of selection bias on instrumental variable (IV) analyses has focused primarily on collider bias. Here, we show that IV estimates can be biased if selection is an effect modifier (or caused by an effect modifier) of the exposure-outcome effect. In a simulation study, we investigate the size and direction of the bias when estimating effects on the mean difference or odds ratio scale. Bias can also occur due to heterogeneity in the instrument-exposure relationship, although its magnitude is typically smaller. Finally, we discuss methods to deal with this bias in practice, including inverse probability weighting (IPW).

The sex ratio at birth (SRB), or secondary sex ratio, is stable, with more boys born alive, but shifts after major events. Acute stress may reduce SRB through disproportionate male fetal loss, while sexual intercourse patterns around conception are associated with offspring sex. We examined whether Brazil's April 2021 Gamma-wave mortality peak was followed by an SRB decline within the previously reported stress-sensitive 3-5 month window. We analysed >81 million live births (1997-2024) from Brazil's Live Birth Information System (SINASC) using pre-pandemic autoregressive forecasts with 95% prediction intervals (PIs). In July 2021, SRB was 50.951% versus 51.185% expected, below the 95% PI (50.968-51.402%), indicating fewer boys. All five macroregions were below expectation, largest in the North, including Amazonas, where Gamma was first identified. Monte Carlo simulation showed that, under usual SRB patterns, a decline with this timing and consistency across Brazil would occur about once in 200 simulations (probability 0.00467), making July 2021 unusual. In February 2023, SRB was 51.441% versus 51.220% expected, above the 95% PI (51.003-51.438%), indicating more boys, nine months after Brazil ended its COVID-19 emergency, suggesting increased sexual intercourse after pandemic relaxation. In April 2023, SRB was 50.956% versus 51.221% expected, below the 95% PI (51.004-51.439%), indicating fewer boys, four months after Pele's death. A comparable UK decline followed Princess Diana's death 4-5 months later. Gamma was prespecified; later observations are hypothesis-generating. These findings highlight SRB's population-level utility and warrant investigation of stress-related and behavioural pathways; ecological design precludes causality.

Accurate prenatal screening from ultrasound imaging is severely constrained by the extreme spatial heterogeneity of fetal structures and the data scarcity inherent to rare congenital anomalies. To overcome these computational bottlenecks, we present RadarMAE, a data-efficient self-supervised representation learning framework tailored for fetal ultrasound analysis. By employing a novel, structure-aware masking strategy, RadarMAE enables robust pre-training on highly limited clinical datasets. Benchmarking the framework on first-trimester screening for trisomy 21 (Down syndrome), our model achieved exceptional diagnostic performance on an internal cohort, with an area under the curve (AUC) of 0.99, 97% accuracy, 94% sensitivity, and 98% specificity. Crucially, through domain-adversarial adaptation, the framework demonstrates robust cross-institutional generalization under extreme source-target data imbalance, bounding performance degradation to less than or equal to 0.04 across all key metrics in external testing. Furthermore, by integrating spatial registration with explainable AI, we constructed the first population-level facial risk atlas for trisomy 21. This data-driven interpretability framework validates standard clinical markers while revealing the comparable diagnostic weight of underutilized features, such as mandibular dysmorphology. Ultimately, this approach provides a scalable, transparent, and high-performance solution for equitable prenatal diagnostics and phenotypic biomarker discovery.

The rapid growth of technologies such as wearable devices and high-throughput sequencing has generated unprecedented volumes of longitudinal data, often with multiple correlated features. This trend highlights the need for robust clustering methods tailored to multi-feature longitudinal datasets. Although previous reviews have benchmarked clustering methods in this context, most evaluations relied heavily on simulated and overly simplified datasets, which substantially limit their interpretability and applicability to real-world data. Missing data is also ubiquitous in longitudinal data, yet understudied in the context of clustering. To address these gaps, we systematically evaluate clustering methods across three distinct types of real-world longitudinal datasets: (1) a high-dimensional, low-timepoint metagenomic study, (2) a low-dimensional, high-timepoint multi-site repeated-measure study, and (3) a high-dimensional, high-timepoint repeated-measure study. We further investigate how missing data and different imputation strategies influence clustering performance. In addition, we introduce a novel resampling framework based on the empirical cumulative distribution function (eCDF) and a copula. This framework enables us to simulate multi-feature longitudinal data that preserves the correlation structure of real datasets, while also allowing for controlled variations such as sample size. Our simulation framework naturally accommodates missing data, thereby providing a flexible tool for benchmarking and comparing clustering methods in more realistic longitudinal settings.

Sample-specific molecular networks, which capture how an individual's gene regulatory structure may deviate from the population as a whole, are important tools for understanding heterogeneity in health and disease. However, most existing methods either do not accommodate uncertainty quantification or cannot handle missing data, non-Gaussian modalities, and longitudinal structure common in modern multi-omic studies. For example, copy number alterations in cancer genomic studies are non-Gaussian and often only partially observed, and their regulatory effects on gene expression may vary substantially across patients, variation that a single, population-averaged network cannot capture. To address these challenges, we propose PRISM (Personalized Relational Inference via Sample-specific Modeling), an empirical Bayes framework for estimating sample-specific covariance matrices from cross-sectional and longitudinal multi-omic data. PRISM models each individual sample's covariance matrix as drawn from an inverse-Wishart prior whose parameters are estimated from the data and allowed to vary across subpopulations or evolve over time, yielding a closed-form posterior and edge-level uncertainty quantification, ensuring scalability to high-dimensional omics data. In simulation experiments, PRISM achieves point-estimate accuracy comparable to or better than existing methods, substantially improved uncertainty calibration under non-Gaussian marginals, and lower computational cost, while maintaining robust coverage under moderate missingness. We apply PRISM to a yeast cell-cycle time course dataset, recovering time-varying biologically meaningful regulatory relationships independently documented in literature, and to paired copy number and gene expression data from ovarian cancer patients, identifying copy-number-expression interactions associated with survival, including TFs HNF4G and PAX8 which have independent evidence of clinical relevance in this cancer.

Spatial transcriptomics enables high-throughput gene expression profiling while preserving spatial information, offering valuable insights into tissue architecture and cellular organization. Methods that achieve single-cell or subcellular spatial resolution typically rely on predefined gene panels, which limits genome-wide discovery. In contrast, sequencing-based spatial transcriptomics platforms provide broad transcriptome coverage but measure expression at lower spatial resolution, capturing mixtures of multiple cell types within each spatial location and thereby requiring accurate cell-type deconvolution. Many existing deconvolution methods do not jointly model molecular expression similarity and spatial context, fail to align predicted cell-type compositions with biological priors, and may produce spatially inconsistent results that do not reflect the smooth and structured organization of real tissues. We present ST-DISTAL, a dual-branch graph convolutional framework for cell-type deconvolution in spatial transcriptomics data. ST-DISTAL integrates conventional spatial graph convolution with spectral Chebyshev filtering through an adaptive attention mechanism, enabling the model to capture both local spatial dependencies and global structural patterns. The framework is trained using a composite objective that jointly enforces accurate cell-type proportion estimation, global distributional alignment with reference profiles, and spatial smoothness. Evaluations on three simulated benchmark datasets demonstrate that ST-DISTAL consistently outperforms state-of-the-art methods in both prediction accuracy and distributional alignment. Further validation on a real human embryonic heart dataset and a colorectal cancer dataset shows biologically plausible and spatially coherent cell-type organization.

Estimating population-level dependence between brain regions is challenging when the number of regions greatly exceeds the number of individuals. We introduce an information-borrowing generalized Bayesian method for high-dimensional covariance estimation from region-level brain measurements. A divide-and-conquer parameterization separates marginal variances from correlation calibrations: an empirical \emph{pilot} correlation matrix--by default, obtained by standardizing a sample covariance matrix (SCM)--provides a reference structure that is adaptively calibrated through region-specific parameters shared across pairwise correlations, reducing the number of free parameters from $\mathcal{O}(d^2)$ to $\mathcal{O}(d)$ ($d$ is the number of regions) while guaranteeing a positive-definite covariance through the inclusion of a ridge term. Horseshoe priors on the calibration parameters enable selective borrowing and shrinkage of pilot information. Under suitable conditions on pilot accuracy, we establish theoretical guarantees for covariance estimation while retaining scalability in ultra-high-dimensional settings via a low-rank likelihood representation and a Woodbury-accelerated Hamiltonian Monte Carlo algorithm. Simulations showed lower estimation error than the SCM, bagging and iteratively refined principal orthogonal complement thresholding (POETRY). Applications to a semi-synthetic cortical-thickness and ADHD resting-state fMRI data revealed significant group differences in the default-mode--visual network and somatomotor network--limbic network blocks after Benjamini--Hochberg correction.